Study of the ionic Peierls-Hubbard model 
using density matrix renormalization group methods 
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Density matrix renormalization group methods are used to investigate the quantum phase dia- 
gram of a one-dimensional half-filled ionic Peierls-Hubbard model at the antiadiabatic limit where 
quantum phonon fluctuations are taken into account partially. We found that two continuous phase 
transitions always exist from dimerized spin-gapped (bond-order-wave) state to band-insulator and 
undimerized spin-gapless (Mott-insulator) phase while undimerized spin-gapless phase vanishes at 
the adiabatic limit. Our results indicate that quantum phonon fluctuations, electron-electron inter- 
action and ionic potential combine in the formation of the bond-order wave phase. 

PACS: 71.30.-fh; 71.10.Fd; 71.10.Pm 



The response of correlated electrons to lattice distor- 
tions in solids has been extensively studied over the years, 
due to its important role in several classes of materials 
including high-Tc cuprates, colossal magnetoresistance 
manganites, conducting polymers and organic charge- 
transfer salts. As a good example, the Peierls-Hubbard 
model, with on-site Coulomb repulsion and lattice dis- 
placement, is a simple yet nontrivial model that exhibits 
a rich ground state phase diagram. Strong correlations 
lead to the separation of charge and spin excitations^ 
while quantum phonon fluctuations can destroy an or- 
dered gapped state^"^. Taking them both into consider- 
ation is essential for a full understanding of the nature of 
these materials. 

With the inclusion of additional terms to the Peierls- 
Hubbard Hamiltonian on different physics background, 
various one-dimensional correlated electronic models 
were actively studied recently such as the ionic Peierls- 
Hubbard model^"^ which is defined as follows 

H = - J2it-(^iui- ui+i)] Bi^i+i + A ^(-l)'n;^ 

].KY,iu,-u,^,f + ±-^Y.P", (1) 

where ni„ is the number operator at site ^, A is electro- 
static potential of cations, and anions in charge-transfer 
salts, ui (pi) is the displacement (momentum) of the site 
Z, a and K are the constant for the electron-phonon cou- 
pling and lattice elasticity, M is the mass and the bond- 
charge density operator -B/,i+i is 
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Here we only considered half- filled case. 

At J7 = 0, model (1) can be solved exactly at the adia- 
batic limit. In Appendix A, we give the detail process and 
discuss the definitions of the quantities we used below. 
Only one phase transition can be found from dimerized 
state which is also called Bond-Order- Wave (BOW) state 
to Band-Insulator (BI) phase as a function of electron- 
phonon (e-p) coupling strength. At A = 0, it is well- 
known that no phase transition will occur even when U 
goes to infinity as long as there exists e-p coupling at the 
adiabatic limit. At U ^ and A ^ 0, earlier work^'^ and 
recent work*"'^ which studied the e-p interaction only in 
the adiabatic limit also concluded that only one phase 
transition, i.e., from the BOW to the BI phase, will be 
present. In other words, the ionic phase with one elec- 
tron per site is always dimerized. Therefore, the phase 
diagram of this model at the adiabatic limit is obvious. 
Only two phases, i.e. BOW and BI phase, can be de- 
tected. 

When the lattice distortion is absent, Eq. (1) repre- 
sents the Ionic Hubbard Model (IHM)^^^^, which was 
used to describe the neutral-ionic phase transition in 
mixed-stack charge transfer crystals^^""'^^ . As pointed 
out by Fabrizio, Gogolin and Nersesyan^, there exists 
an unusual spontaneously dimerized insulator phase, the 
BOW phase, which separates the BI from the Mott in- 
sulator (MI) phase^. However, in reality, a lattice dis- 
tortion always exists and it couples to electronic degrees 
of freedom strongly in these crystals. Structure changes, 
such as volume contraction could be used as an external 
parameter to drive the neutral-ionic transition^*. Pho- 
toinduced cooperative phenomena were also observed^^. 
So an important issue to address is the effect of electron- 
phonon interactions on the phase diagram. 
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However, it is well known that results obtained at 
the adiabatic limit are unreliable. For the Su-SchriefFer- 
Heeger (SSH) model, Fradkin and Hirsch^'^ pointed out 
that the low-energy behavior of the system is actually 
governed by the antiadiabatic limit M = 0, rather than 
the adiabatic limit M oo. The system at any non- 
zero frequency is renormalized to the limit of infinite fre- 
quency. For the Holstein modeP and the spin-Peierls 
model'', more sophisticated calculations showed that uni- 
form gapless phase exists unless the e-p coupling is suf- 
ficiently large. On the other hand, the system is found 
always in the dimerized gapped state at the adiabatic 
limit. Therefore, in order to study the effect of e-p inter- 
action truly and understand the whole phase diagram, it 
is necessary to investigate the ground state properties of 
the system at the antiadiabatic limit. 

In the present work, we perform an extensive numerical 
study of Eq. (1) at the antiadiabatic limit using the den- 
sity matrix renormalization group (DMRG)^^ technique. 
We found that, even with very strong e-p coupling, two 
continuous phase transitions from BI phase to MI phase 
are obtained from this model. One is the spin transi- 
tion ai U = Us where spin excitation gap closed, and the 
other is the charge transition at U = Uc < Us where the 
charge excitation gap vanishes. Between these two crit- 
ical points, the system is dimerized. In contrast to the 
adiabatic limit where MI phase vanishes^"^, this means 
the ionic phase can also be undimerized after considering 
quantum phonon fluctuations. Our phase diagram is sim- 
ilar to that of IHM^^'^^ which is still highly controversial 
over the years'^^^'"'. Furthermore, in the region of strong 
ionic potential, the critical values Uc and Ug decrease si- 
multaneously with increasing e-p coupling, while Us will 
increase with increasing e-p coupling at sufficiently small 
ionic potential. The ground state phase diagram is ob- 
tained with the use of finite-size-scaling analysis. 

At the antiadiabatic limit M = 0, an effective inter- 
acting fermion model^'' 

H = -tY,Bi,i+i + UY^ (ni^ - I) («U - \) 

-W^{B^+,f+A Y^i-lYni,, (3) 

; /cr 

can be obtained where the effective bond-charge attrac- 
tion {W = /2K) term accounts for the contribution of 
the phonon quantum fluctuations. So the Hamiltonian 
(3) could be viewed as a ID e-p interacting system in- 
cluding both the quantum phonon fluctuations and the 
electron correlations. 

In this paper, we have applied the finite-size DMRG 
algorithm with open boundary conditions to study the 
Hamiltonian (3) at half-filling. This method allows us 
to probe directly correlation functions and structure fac- 
tors associated with the spin density wave (SDW), the 
charge density wave (CDW) and the bond order wave in 
the ground state. Lattices up to 512 sites were frequently 
used in our studies. The largest number of states kept 



in the calculation was m = 512 per block. The hopping 
integral t is set to 1 as the energy unit. The weight of 
the discarded states was typically about 10~^ — 10""'^'^ de- 
pending on whether the system is in its critical state or 
not in the final sweep. The convergence tests as functions 
of number of states kept were carefully performed. We 
checked our DMRG calculations against exact numerical 
results for noninteracting {U = W = 0) chains (up to 512 
sites) and results from exact diagonalization for interact- 
ing (U 0, W 0) chains (up to 14 sites). Excellent 
agreement was found in both cases. When interactions 
are turned on, there exist finite excitation gaps on finite 
chains, so the accuracies of all quantities we calculated 
are no worse than that of the noninteracting case. Thus, 
numerical errors in our work could be safely estimated to 
be smaller than 10~^. 




FIG. 1. The staggered BOW correlation functions of the 
ionic Peierls-Hubbard model for A/t = 0.30 and W/t = 0.30. 

The inset shows an extrapolation of m%ow {^) with a 
third-order polynomial in 1/L. 

Now let's introduce a way to determine the phase 
boundaries accurately. Although parts of this method 
have already been first used in determining the phase 
boundaries of IHM^^, we perform this method more care- 
fully and systematically in this paper than ever by using 
various finite-size analysis. First we defined the stag- 
gered BOW correlation function which is the most direct 
evidence for the long-range BOW state 

CbOW (r) = (-l)"" (y^ Yl {Bl,l+lBl+r,l+r+l) - e") 

(4) 

where B = Ylii In Fig. 1, we show the stag- 

gered BOW correlation functions with increasing Hub- 
bard C/ at A/t = 0.30 and W/t = 0.30. Here we only 
perform the average in eq. (4) over 256 sites in the mid- 
dle of the 512-site system, i.e. the sum over I is from 129 
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to 256 and Lav = 128, in order to further avoid bound- 
ary effects. L is the half of the chain length N which 
is located in the central of the system to further elimi- 
nate edge effects. The results indicate that there exist 
three different phases in model (3) since the staggered 
BOW correlation functions show three distinct type of 
behavior as r increases: (i) it decays exponentially at 
U/t = 0.50, indicating that the system has no BOW or- 
der; (ii) it converges to a nonzero constant at U/t = 1.45 
and U/t = 1.55, indicating that the system is in the 
BOW phase with a finite width; (iii) it decays as 1/r 
at U/t = 2.50, indicating that the system is in another 
phase. The BOW order parameter in the thermodynamic 
limit 



(6) 



^BOW 



lim 

L—^oo L 



(5) 



can be obtained by fitting 



AL) i=iErC'BOw{r), 



BOW 

where the sum over r goes to L) with a third-order poly- 
nomial in 1/L since m^Q^ (L) ^%ow fo'" L ^ oo. 
The inset of Fig. 1 shows such extrapolations. In the 
following we always take Lav — 4 which is enough to 
minimize the oscillations due to the open boundary con- 
ditions and sum over r is only take over central L site 
where L is half of the chain length N to further elim- 
inate edge effects. We find that m']^ow i-^) approac;lK>s 
zero when U/t = 0.50 and 2.50 but remains finite when 
U/t = 1.45 and U/t = 1.55 which indicate that there 
exists a finite region of the BOW phase. 




FIG. 2. Behavior of S'sow(7r) across the phase boundary 
for A/t = 0.3, W/t = 0.30. The inset shows a linear extrap- 
olation of the critical values Uc and Us with 1/L. 

Next we have studied the nature of BI-BOW and 
BOW-MI transitions by calculating the static structure 
factors corresponding to different phases to determine the 
phase boundaries. The first structure factor studied is 
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According to Fabrizio et. al.®, phase transitions on the 
BI-BOW and BOW-MI phase boundaries are an Ising 
type and KT type respectively, the staggered connected 
correlation function falls off algebraically as 



( — 1)'^ {{Bi^l^lBi+r,l+r+l) — (BlJ+l) {Bi+r,l+r+l)) 



(7) 
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FIG. 3. Finite-size analysis for 5soM'(7r) and ScDw{Tr) in 
the vicinity of the first phase transition at Uc for A/t = 0.30, 
W/t = 0.30. 

Away from phase boundaries, this quantity falls off 
exponentially. Therefore the Sbow M is expected to 
diverge at these two critical points if ry < 1 or reach max- 
imums if > 1 as the system size goes to infinity. Fig. 
2 shows the results of the Sbow (ti") for different system 
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sizes with A/t = 0.30, W/t = 0.30 and {)<U/t< 3. As 
expected, tlie Sbow {t^) peaks twice for all the different 
system sizes we calculated. The positions of these two 
peaks become closer as the system size is larger. The 
inset of Fig. 2 shows linear extrapolations of the posi- 
tions of these two peaks with l/L. We find that these 
two peaks will not merge at L — » oo which indicate again 
that the BOW phase remains finite at thermo dynamic 
limit. In order to give more convincing evidence, we did 
another finite-size analysis in the vicinity of these two 
phase transitions. Let us start from the first phase transi- 
tion &iU = Uc- Fig. 3(a) presents plots of \n[SBOw ('"■)] 
versus ln[L] for A/t = 0.30, W/t = 0.30 and three dif- 
ferent values of U /t around the first critical point. Data 
points for U/t = 1.26 indeed fall on a straight line, indi- 
cating critical scaling for the BOW fluctuations. At the 
other two points U/t — 1.20 and 1.30, data points bcliave 
nonlinearly due to the exponential decay term. Applying 
the same finite-size analysis to the CDW structure fac- 
tor ScDW (^)j wc can also explore the nature of the first 
phase transition. The CDW structure factor is defined 
as 



SsDW (?) = T — 



^CDW 



-Lj2e^'^^{{nm+r)-{ni){ni+r)) (8) 



The linear behavior of ln[ScDW (tt)] aroimd U/t = 1.26, 
shown in Fig. 3(b), confirms the vanishing of the charge 
gap at the first phase transition point. 



Ir 



(9) 



It is well known^^ that if the ground state of a ID sys- 
tem is spin-gapless, the spin-spin correlation falls alge- 
braically with exponent equal to 1. It has been further 
shown^^ that in the spin-gapless phase Ssdw (<z) /q 
I/tt as (7 — !■ whereas in the spin-gapped phase 
Ssdw {q) /q 0. Even a very small spin gap can be 
detected in this way, since it is in practice sufficient to 
see the ttSsdw (s) /q decay below 1 for small q to con- 
clude that a spin gap must be present. Fig. 4 shows the 
behavior of ttSsdw {q) /q for A/t = 0.30, W/t = 0.30 
and different values of U/t. In the gapless region, loga- 
rithmic corrections^^ make it difficult to observe the ap- 
proach to 1 as — ^ 0. In analogy with spin systems^®, 
we expect the leading logarithmic corrections to vanish 
at the point where spin gap opens and therefore exactly 
at the critical point there should be a clear scaling to 1. 
Based on results shown in Fig. 4, we estimate the MI- 
BOW boundary to be at U/t = 1.70±0.02 at A/t = 0.30, 
W/t = 0.3 which is consistent with the results shown in 
Fig. 2. The inset of Fig. 4 further provides evidence for 
the transition from the spin-gapped state, as identified 
by the exponential decay of the staggered SDW correla- 
tion function, to the spin-gapless state, as characterized 
by the 1/r-decay of the SDW correlation function. 



1,0- 



0,5 



0,0 




U/t=0.50, N 


=256 


10' 


U/t=1.45, N 


=256 


U/t=1.60, N 


=256 


1 


U/t=1.70, N 


=256 


U/t=1.80, N 


=256 


10 


U/t=2.50, N 


=256 




U/t=0.50, N 


=512 


10- 


U/t=2.50, N 


=512 





U/t=1.45 
U/t=2.50 
-1/r fitting 




0,0 



1,2 



2,0 



2,8 



FIG. 4. ■KSsDw{q)/q vs q for A/t = 0.30, W/t = 0.30 and 
6 different values of U across the MI-BOW boundary. 

Next wc determine the nature of the second phase tran- 
sition at U = Us- As predicted by Fabrizio et. al.^, it is a 
quantum phase transition of the KT type. This makes it 
difficult to determine the phase boundary directly from 
the behavior of Sbow (tt) and Ssdw ("■) (defined below) 
due to the finite-size effects. Instead, we apply an in- 
direct method, used by Scngupta et al.**, to confirm the 
second phase transition. The SDW structure factor is 
defined as 



CsDW {r) = j{-irY.('i'Ur) 



(10) 
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FIG. 5. Phase diagram of model (3) for two different e-p 

couplings. The inset shows different behavior of Us for weak 
(A/t = 0.10) and strong (A/t = 0.30) ionic potentials as 
a function of e-p coupling W. While the beiiavior of Uc as 
a function of W is always the same for all ionic potentials. 
Here we only show the plot for Uc at A/t = 0.30. 
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Finally, wc present in Fig. 5 the resulting phase di- 
agram in the [/ — A plane for two values of W. For 
A = U = 0, model (3) becomes the t — W model^'' 
which can be mapped from the SSH model at the antiadi- 
abatic limit. This model has been studied by the DMRG 
method^^ and the renormalization group analysis^^. The 
ground state is dimerized and the BOW order parameter 
is nonvanishing as long asW^O. After switching on the 
Hubbard [/ or a finite ionic potential A, the BOW phase 
could be destroyed. The model undergoes quantum phase 
transitions from the BOW phase either to the MI phase 
or the BI phase. The critical value U and A will increase 
with increasing e-p coupling W. In the weak ionic po- 
tential region, such as A/t = 0.1, on increasing the e-p 
coupling W, the transition points Uc and Us move apart 
and the separation between Uc and Ug becomes signifi- 
cantly larger. However, for strong ionic potential, such 
as A/t = 1.0, both Uc and Us decrease and the width 
of the BOW phase increases slightly with increasing e- 
p coupling. When the ionic potential is intermediate, 
such as A/t = 0.3, Uc and Us will also decrease simul- 
taneously while the width of the BOW phase increases 
significantly. For model (1) at the adiabatic limit, a spin 
gap will always be present*. However, at the antiadia- 
batic limit, the transition from BOW to MI phase always 
occurs even though the e-p coupling is sufficiently large. 

The inset of Fig. 5 shows the asymptotic behavior 
of the critical points Uc and Us at A/t = 0.10 and 



A/t = 0.30 with increasing c-p coupling. For weak on- 
site potential such as A/t = 0.10, the critical point Us 
becomes larger with increasing c-p coupling while for 
strong on-site potential for example A/t = 0.30, it de- 
crease monotonously with the increasing e-p coupling. 
However the behavior of Uc is the same for all value of 
on-site potential, it decreases monotonously as a func- 
tion of c-p coupling. Here we only show one set of data 
at A/t = 0.30. 

In conclusion, we have studied ID half-filled Ionic 
Peierls-Hubbard model at the antiadiabatic limit using 
the DMRG method. The phase diagram is obtained by 
investigating correlation functions and structure factors. 
In contrast to the adiabatic limit, the transition from the 
dimerized spin-gapped state to the spin-gapless state oc- 
curs for any value of the ionic potential A. Compared to 
the IHM, the BOW phase always exists in the presence of 
e-p coupling W. The first critical value Uc for the charge 
gap always decreases with increasing W while the second 
critical value Us for the spin gap shows different behavior 
for weak and strong ionic potential. 
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APPENDIX A: EXACT SOLUTION TO MODEL (1) AT THE ADIABATIC LIMIT IN THE 

NONINTERACTING CASE 



In this appendix we solve model (1) exactly at the adiabatic limit in the noninteracting case. First we define the 
dimensionless coupling constant as A = 2a'^/Kt and dimerized order parameter as S = au/t where (—1) u = u; — u;+i. 
Then model (1) can be rewritten as following 



H/t = -j2[i- (-1)' s] Bi,i+, + D j2{-iym^ + jT.^' 



(Al) 



where D = A/t. Using Fourier transformation and unitary transformation, the original operator can be expressed 



as 



^1 = 4^12^''' Hi04.+(3tii)bt) 



(A2) 



where 



and k G (-|, |]. Here 



where 



oik (0 =Uk + (-1)' Vk, Pk (0 = -^fe + (-1)' Wfe 



\/Ek - 2 cos A; 
"fe = 7^W= ' ^fe 



—2i6smk 



^2Ek {Ek- 2 cos k) 



(A3) 



(A4) 



Ek = V4cos2fc-F452 sin^ k + D^. 



(A5) 
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Then the Hamiltonian (Al) can be diagonahzed as 

H/t = Y.Ek [ataua - b+^bk.) . (A6) 

The quantities related to the BOW phase can be expressed as following 

{Bi+r,i+r+i) = 4 (Ai - (-!)'+'■ Bi) . (A7) 
here Bi is the order parameter, Ai is the average bond length and 

{Bl+r,l+r+lBlJ+l) = {Bi+r,l+r+l) (Bl^l+l) + 



- (Ar+i - {-if Br+l) + (-1)' Br-l) 

4[Cr+iCr-i-{A^,-B^^)],r = odd 



, r = even 



(A8) 



where 



1 /■ 2 „ 2 cos kr cos k , . „, 

- / dk . (A9) 



„ S „2smkrsmk . 
Br = - dk . (AlO) 

7!" Jo Ek 

^ D COSk /A-,-,^ 

Cr = — dk——. (All) 

TT Jo Ek 

Then we find that the definition (4), (5), and (6) in our paper is usefuU. Since B = lT,i{Bi,i+i) = 4Ai, we can 
obtain 



CBOwir) = 16B? + 

4 [C^ - {Ar+iAr-i + Br+iBr-i)] , T = eveu 
4 - B^) - Cr+iCr-i] ,r = odd 



(A12) 



Obviously, if the system is in the BOW phase, CBOwir) 16Bf remains constant at large r otherwise Csowir) — * 
in the CDW phase since second line in eq. (A12) is exponential decay with the distance r in both phases. The BOW 
order parameter defined in our paper can be obtained as 

Abow = lim y V (-1)'+' = 4Bi. (A13) 

L— >oo Li ^-j'' 

Meanwhile 

m%ow (oo) = lim y V CBOw{r) = 16Bf (A14) 

r 

since the second line in eq. (A12) is exponential decay with the distance r. Here we should mention that on the phase 
boundary the second line in eq. (A12) is power-law decay. Nevertheless eq. (A14) is still valid. The staggered BOW 
structure factor is 

Q (-^\ ( 'i[C^ - iAr+iAr-i+ Br+iBr-i)] ,r = even\ 

SBOw{n) = ^ ' ^[{A^^_B^^)-Cr+iCr-i],r = odd j' ^^^^^ 

As we mentioned above, exactly on the phase boundary it is power-law decay otherwise; exponential decay. Finally 
we give the self-consistent equation which can be used to determine the phase diagram in the A — A plane 

2A /"2 „ 2sinfcrsinfc z.^^n 
dk — . (-^16) 

Ek 
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